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Abstract 

Calculations of observables in Quantum Chromodynamics are typically per- 
formed using a method that combines numerical integrations over the mo- 
menta of final state particles with analytical integrations over the momenta 
of virtual particles. I discuss a method for performing all of the integrations 
numerically. 
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This Letter concerns a method for performing perturbative calculations in Quantum 
Chromodynamics (QCD) and other quantum field theories. Specifically, I am concerned 
with cross sections and other QCD observables in which one measures something about the 
hadronic final state. Here one cannot use the special techniques that apply to inclusive 
quantities like the structure functions in deeply inelastic lepton scattering. The general 
class of calculations of interest in this Letter includes jet cross sections in hadron-hadron and 
lepton-hadron scattering and in e + e~ — * hadrons. There have been many calculations of this 
kind carried out at next-to-leading order in perturbation theory, resulting in an impressive 
confrontation between theory and experiment and in an accumulation of evidence supporting 
QCD as the correct theory of the strong interactions J^j. These calculations are based on 
a method introduced by Ellis, Ross, and Terrano || in the context of e + e~ — > hadrons. 
Stated in the simplest terms, the Ellis-Ross- Terrano method is to do some integrations over 
momenta Pi analytically, others numerically. I shall argue that it is possible instead to do 
all of these integrations numerically. Furthermore, I shall argue that performing all of the 
integrations numerically has some advantages, principally in the flexibility that it allows. 

In this Letter, I address only the process e + e~ — > hadrons. Thus I do not address 
the issue of factorization that is associated with initial state hadrons. I discuss three-jet- 
like infrared safe observables at next-to- leading order, that is order a 2 . Examples of such 
observables include the thrust distribution, the fraction of events that have three jets, and 
the energy-energy correlation function. 

Let us begin with a precise statement of the problem. The order a 2 contribution to the 
observable being calculated has the form 



Here the do^} are the order a 2 contributions to the parton level cross section, calculated 
with zero quark masses. Each contains momentum and energy conserving delta functions. 
The da^ include ultraviolet renormalization in the MS scheme. The functions S describe 
the measurable quantity to be calculated. We wish to calculate a "three-jet" quantity. 
That is, S2 = 0. The normalization is such that S n = 1 for n = 2,3,4 would give the 
order a 2 s perturbative contribution the the total cross section. There are, of course, infrared 
divergences associated with Eq. ([]]). For now, we may simply suppose that an infrared cutoff 
has been supplied. 

The measurement, as specified by the functions S n , is to be infrared safe, as described 
in Ref. ||: the S n are smooth functions of the parton momenta and 



for < A < 1. That is, collinear splittings and soft particles do not affect the measurement. 




(1) 



S n+ i(k h . . . , Xk n , (1 - \)k n ) = «S n (fei, ...,k n ) 



(2) 
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It is convenient to calculate a quantity that is dimensionless. Let the functions S n be 
dimensionless and eliminate the remaining dimensionality in the problem by dividing by a , 
the total e + e~ cross section at the Born level. Let us also remove the factor of (o> s /tt) 2 . 
Thus, we calculate 

a® 

1= - ( i , 2 - (3) 

We note that X is a function of the cm. energy yfs and the MS renormalization scale \x. 
We will choose fi to be proportional to yfs: \i = A UV y/s. Then X depends on A uv . But, 
because it is dimensionless, it is independent of yfs. This allows us to write 



X = / dyfi h(y/s~) I{A UV , yfi), (4) 
J 

where h is any function with 

/ d^~s h(y/s) = 1. (5) 

J 

The integration over yfs eliminates the energy conserving delta function in X. The physical 
meaning is that, by smearing in the energy yfs, we force the time variables in the two current 
operators that create the hadronic state to be within 1 / y/s of each other. Thus we have a 
truly short distance problem. 

I now describe how one would calculate X using the Ellis-Ross- Terrano method. Each 
partonic cross section in Eq. (|TJ) can be expressed as an amplitude times a complex conjugate 
amplitude. One must calculate the amplitudes in 4 — 2e dimensions. (In the case of the 
process e + e~ — > hadrons, this calculation was performed by Ellis, Ross, and Terrano ||.) 
For tree diagrams, the calculation is straightforward. For loop diagrams, this involves an 
integration, which is performed analytically. The integrals are divergent in four dimensions, 
so one obtains divergent terms proportional to 1/e 2 and 1/e in addition to terms that are 
finite as e — > 0. Having the amplitudes and complex conjugate amplitudes, one must now 
multiply by the functions S n and integrate over the final state parton momenta. These 
integrations are too complicated to perform analytically, so one must use numerical methods. 
Unfortunately, the integrals are divergent at e = 0. Thus one must split the integrals into 
two parts. One part can be divergent at e = but must be simple enough to calculate 
analytically. The other part can be complicated, but must be convergent at e = 0. One 
calculates the simple, divergent part and cancels the 1/e 2 and 1/e pole terms against the 
pole terms coming from the virtual loop diagrams. This leaves the complicated, convergent 
integration to be performed numerically. 

This method is a little bit cumbersome, but it works and has been enormously successful. 
However it has proven to be difficult to apply the method in the case of two virtual loops. 
Even with one virtual loop, the method is not very flexible. Any modification of the integrand 
requires one to recalculate the amplitudes, and the modification must be simple enough that 
one can calculate the amplitudes in closed form. 

Let us, therefore, inquire whether there is any other way that one might perform such a 
calculation. We note that the quantity X can be expressed in terms of cut Feynman diagrams, 
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FIG. 1. Two cuts of one of the Feynman diagrams that contribute to e + e — * hadrons. 

as in Fig. |I[ The part of the diagram to the left of the cut is a term in the amplitude. The 
part to the right of the cut is a term in the complex conjugate amplitude. The dots where 
the parton lines cross the cut are intended to represent the function S n (k\, . . . , k n ). Each 
diagram is a three loop diagram, so we have integrations over loop momenta £^ and £%. 

We first perform the energy integrations. For the graphs in which four parton lines cross 
the cut, there are four mass-shell delta functions 5(kj). These delta functions eliminate the 
three energy integrals over £®, £\, and £\ as well as the integral (PS) over sfs. For the graphs in 
which three parton lines cross the cut, we can eliminate the integration over yfs and two of 
the £°j integrals. One integral over the energy E in the virtual loop remains. The integrand 
contains a product of factors 



E - Q°j - coj + xe E - Q°j + oo,j - ie w 

where E — Q°j is the energy carried by the Jth propagator around the loop and ujj is the 
absolute value of the momentum carried on that propagator. We perform the integration 
by closing the integration contour in the lower half E plane. This leads to n terms for a 
virtual n point subgraph. In the Jth term, the propagator energy E — Q°j is set equal to the 
corresponding uj and there is a factor l/(2uj). Note that the entire process of performing 
the energy integrals amounts to some simple algebraic substitutions. 

Let us denote the contribution to X from the cut C of graph G by X{G,C). This 
contribution has the form 

X(G, C)= J d% d% d 3 £ 3 g{G, C; £), (7) 

where £ denotes the set of loop momenta {£\, £2^3}- The functions g have some singularities, 
called pinch singularities, that cannot be avoided by deforming the integration contour and 
some non-pinch singularities that can be avoided by a contour deformation. A detailed 
analysis, which I sketch below, has been given by Sterman 

The pinch singularities occur when one parton branches into two partons with collinear 
momenta or when one parton momentum goes to zero. These singularities can lead to 
logarithmic divergences in the corresponding integral. (As pointed out in the q^q u terms 
in gluon self-energy subgraphs lead to quadratic divergences. I will discuss these terms 
shortly.) If we were to calculate the total cross section by using measurement functions 
«S n = 1 in g, then the singularities would cancel M between the functions g associated 
with the various cuts C of the same graph G. The underlying reason is unitarity. Now 
in our case of measured shape variables, the values of S n corresponding to different cuts 



4 



C are different. This would ruin the cancellation, except that just at the collinear or soft 
points the functions S n match. Thus the singularities present in the individual g(G,C;£) 
cancel in the sum, 2c9(G, C\ £). To be precise, the cancellation reduces the strength of the 
singularity from a strength sufficient to give a logarithmically divergent integral to one that 
gives a convergent integral. 

The function g also has singularities that can be avoided by deforming the integra- 
tion contours, the non-pinch singularities. These can be characterized rather simply in the 
present case of a single virtual loop. Let £ be the a loop momentum that flows around this 
virtual loop. We may choose £ so that it is the momentum carried by one of the loop prop- 

— * 

agators and so that a momentum Q flows out of the loop and into the final state between 
the propagator in question and another propagator further along the loop. This second 
propagator carries momentum £ — Q. A positive energy Q° also flows out of the loop and 
into the final state. Since (Q°, Q) is the four-momentum of a group of final state parti- 
cles, we have Q° > \Q\. Then there is a singularity when \£\ + \£ — Q\ = Q°. In the case 
Q° = \Q\, this is a collinear singularity, and will cancel between cuts. In the case Q° > \Q\, 
the singular surface is an ellipsoid. The singularity does not cancel, but the Feynman rules 
provide an it prescription that tells us that we should deform the £ integration contour 
into the complex £ plane so as to avoid the singularity. Here deforming the contour means 
replacing £ by a complex vector £ + m. Then one simply chooses the imaginary part, k, of 
the loop momentum as a function of the real part, £, and supplies an appropriate jacobian 

— » 

J . Since the momentum £ that flows around the virtual loop in question is, in general, a 
linear combination of the three loop momenta £j, one should write the general relation as 

— * — # 

£j — > £j + iKj, or, in a abbreviated notation, £ — > £ + in. Then 

X(G, C) = j d% d% d% J{G, C; £) g{G, C;£ + m(G, C;£)). (8) 

Note that the contour deformation is more than just a mathematical trick. The analyticity 
that allows it is a consequence of the causality property of the field theory. 

Now, there is a simple argument that the pinch singularities cancel between different cuts 
C of the same graph G. There is another simple argument that the non-pinch singularities 
are not a problem because one can deform the integration contours so as to avoid them. If 
one is being careful in the proof, one notes that the deformations required to escape the non- 
pinch singularities are different for each cut C . (The deformations for different cuts must be 
different if one wants the momenta of final state particles to be real.) On the other hand, 
the cancellation of pinch singularities requires that the contours for different cuts C be the 

same as one approaches the pinch singularities. That is, the set of vectors {ki, k 2 , k 3 } that 

— » — * — » 

define the deformation associated with a cut C must to go to zero as {£±,£2, £3} approaches 
a pinch singular surface. 

This sort of consideration of how to prove that X is finite instructs us how to actually 
calculate X. With the contours appropriately chosen, the integral 

1(G) = f d% d% d% ]T J(G, C; £) g(G, C;£ + m(G, C; £)) (9) 
J c 

is finite. One can simply compute it by Monte Carlo integration. Note the significance of 
putting the summation over cuts inside the integral. When we sum over cuts for a given 
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point in the space of loop momenta, the soft and collinear divergences cancel because the 
cancellation is built into the Feynman rules. If we were to put the sum over cuts outside 
the integration, as in the Ellis-Ross- Terrano method, then the individual integrals would be 
divergent. The calculation would thereby be rendered more difficult. 

Eq. (^) represents the main point of this Letter. It is perhaps useful to point out that one 
must remove from the integration tiny regions near the collinear and soft singular points, 
since otherwise roundoff errors would spoil the cancellation of the individual contributions. 
I leave for later, more specialized, papers the details of how one can choose a specific contour 
deformation so that the theoretical cancellation is realized in practice. I also skip a discussion 
of how one can choose the density of integration points so as to perform the numerical 
integration by the Monte Carlo method. There are, however, two points of principle that 
I discuss below, since they are not analyzed in Ref. M: renormalization and the special 
treatment required for self-energy subgraphs. 

Renormalization is conventionally accomplished by performing loop integrations in 4 — 2e 
space-time dimensions and subtracting the resulting pole term, c/e. Clearly that is not 
appropriate in a numerical integration. However, one can subtract instead an integral in 
4 space-time dimensions such that, in the region of large loop momenta, the integrand of 
the subtraction term matches the integrand of the subgraph in question. The integrand of 
the subtraction term can depend on a mass parameter \x in such a way that the subtraction 
term has no infrared singularities. Then, one can easily arrange that this ad hoc subtraction 
has exactly the same effect as MS subtraction with scale parameter \i. 

Self-energy subgraphs require a special treatment. Consider a quark self-energy sub- 
graph — zS with one adjoining virtual propagator and one adjoining cut propagator. This 
combination really represents a field strength renormalization for the quark field, and is 
interpreted as 



q z 



2n5(q 2 ). (10) 

q 2 =0 



In order to take the q 2 — > limit here while at the same time maintaining the cancellation 
of collinear divergences, we write 



^{q)4 g 



2 



q 2 (2tt) 



C F / d 3 £ 
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|)2_ (g0) 2- 



(11) 



— » * _ -* 

where k± = ^q± I. This expression is obtained by writing the left-hand side as a dispersive 
integral with the cut self-energy graph appearing as the discontinuity. When the virtual self- 
energy is written in this form, the g 2 — *> limit is smooth and, in addition, the cancellation 
between real and virtual graphs in the collinear limit £ cx qis manifest. It should be noted 
that the integral in Eq. (|TT|) is ultraviolet divergent and requires renormalization, which can 
be performed with an ad hoc subtraction as described above. 

One may expect that the representation of virtual propagator corrections in terms of the 
cut propagator will prove to be convenient in future modifications of the method described 
in this Letter. One may want to make modifications to the gluon propagator, in particular, 
in order to implement a running strong coupling and to insert models for the long distance 
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propagation of gluons. In addition, one may want to modify propagators so that partons 
can enter the final state with q 2 > so as to mesh an order a 2 perturbative calculation with 
a parton shower Monte Carlo program. 

Consider, finally, the one loop gluon self-energy subgraph, 7r M1/ (g). The term in ir fll/ (q) 
proportional to q^q u contributes quadratic infrared divergences [Q. This problem can be 
addressed by replacing by P^7i al3 Pg, where Pg = g£ — q^q a /q 2 , with q = (0,<f). The 
terms added to it^ are proportional to either g M or q u and thus vanish when one sums over 
different ways of inserting the dressed gluon propagator into the remaining subgraph. Since 
P£q a = 0, the problematic q^q u term is eliminated. Effectively, this is a change of gauge for 
dressed gluon propagators from Feynman gauge to Coulomb gauge. 

We have seen that a completely numerical integration of the cut Feynman diagrams for 
a physical quantity can, in principle, produce the numerical value of the quantity. Further- 
more, there may be advantages in simplicity and flexibility associated with this approach. 
The question naturally arises, can such a calculation work, in a practical sense? In order 
to find out, I have constructed a demonstration computer program along the lines out- 
lined above 0. I have used the program to calculate da^/dT, the order a 2 contribution 
to the thrust distribution. More precisely, I have calculated the ratio R(T) of da^/dT to 
{da^/dT} K N, where [da^/dT] KN is a fit to the tabulated results for da^/dT as given by 
Kunszt and Nason 0. In the range 0.71 < T < 0.95, the function [da^/dT]^ varies by 
about a factor of 30. The ratio R(T) should be 1. The results are reported in Fig. |^. We 
can conclude that the automatic cancellations between different cuts of the same diagram 
are indeed realized and that completely numerical integration for QCD observables beyond 
the leading order is a practical possibility. 

Outlook. Substantial effort will be required to test the computer code used for Fig. || in 
order to detect and remove any bugs that it may contain and then to document the code 
and the algorithms used. This work will be reported in future papers. With this code in 
hand, or with improved code from other authors, one can attack more difficult problems than 
discussed here. It remains to be seen for what problems the completely numerical method 
will prove to be more powerful than the analytical/numerical method that has served us so 
well up to now. 

I thank U. Amaldi for initial encouragement to somehow do better with QCD calcula- 
tions, L. Surguladze for advice on Feynman numerators, Z. Kunszt and P. Nason for help 
with the comparison to the results ||, and T. Sjostrand for advice about future applications. 
This work was supported by the U. S. Department of Energy. 
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FIG. 2. Ratio R of the order a 2 s contribution to the thrust distribution as calculated here to 
the same quantity as calculated by Kunszt and Nason, Ref. ||. The horizontal lines represent 
the expected result, 1, with an error estimate. The error bars on the computed points represent 
statistical errors. There is some correlation expected between neighboring points. 
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